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Summary 

The full-approximation-scheme (FAS) multigrid methc d is applied to several implicit flux-split algorithms 
for solving the three-dimensional Euler equations in a b<idy-fitted coordinate system. Each of the splitting 
algorithms uses a variation of approximate factorization iind is implemented in a finite-volume formulation. 
The algorithms are all vectorizable with little or no scab r computation required. The flux vectors are split 
into upwind components using both the Steger- Warming and Van Leer splittings. Results comparing pressure 
distributions with experimental data using both splitting rypes are shown. The stability and smoothing rates 
of each of the schemes are examined using a Fourier analysis of the complete system of equations. Results 
are presented for three-dimensional subsonic, transonic, a nd supersonic flows that demonstrate substantially 
improved convergence rates with the multigrid algorithm The influence of using both V-cycle and W-cycle 
strategies on the convergence is examined. By using the multigrid method in both subsonic and transonic wing 
calculations, the final lift coefficient is obtained to within (!.l percent of its final value in as few as 15 multigrid 
cycles for a mesh with over 2 10 000 points. A spectral radiu s of 0.890 is achieved for both subsonic and transonic 
flows over the ONERA M6 wing, whereas a spectral radius of 0.830 is obtained for supersonic flow over an 
analytically defined forebody. Results compared with exp iriment show good agreement for all cases. 

Introduction 

Upwind difference methods for solving the Euler equa tions are becoming increasingly popular for several 
reasons. The time-dependent Euler equations form a systi m of hyperbolic equations, and upwind differencing 
models the characteristic nature of the equations in that information at each grid point is obtained from 
directions dictated by characteristic theory. Upwind methods have the advantage of being naturally dissipative. 
Separate spatial dissipation terms, such as those generally required in a central difference method to overcome 
oscillations or instabilities arising in regions of strongly vr rying gradients, need not be added. Some examples 
of upwind methods include the A-method (ref. 1), the spli -coefficient method (ref. 2), the flux-vector-splitting 
method (refs. 3, 4, and 5), and the flux-difference-splittin,! ; method (ref. 6). 

Although the A-method and the split-coefficient methc.d closely mimic the method of characteristics, they 
are both applied to the nonconservative form of the equations and consequently require the use of shock-fitting 
techniques to obtain the correct location and strength of sl ocks in transition flows. Use of the conservationdaw 
form of the Euler equations allows shock waves to be caj>tured as weak solutions to the governing equations 
and circumvents the difficulty in applying shock-fitting te. hniques to arbitrary flows. Both the flux-difference- 
splitting and flux-vector-splitting methods can be applied to the conservation-law form. 

The upwind method used in the current work is the flux-vector-splitting method in which the flux vectors 
are split into forward and backward contributions based on an eigenvalue decomposition and are differenced 
accordingly. The splittings investigated include those o Steger- Warming (refs. 3, 4, and 7) and Van Leer 
(refs. 5, 8, and 9). In comparison with unsplit methods, the advantages of flux splitting are obtained at the 
cost of increased computational work since two sets of fluxes are computed for each coordinate direction and 
implicit schemes require two sets of flux Jacobians for cciisistent linearization of the fluxes. In addition, the 
split fluxes and flux Jacobians are also generally more cor iplicated than the unsplit terms because of the logic 
involved with eigenvalue sign changes. 

In order to offset the additional computational worl; of the upwind methods, it is highly desirable to 
accelerate the convergence rate, especially when only steady-state solutions are sought, the objective is to 
reduce the computer time required while still maintaining the high level of robustness and accuracy attained 
from upwind differencing. Accelerating the convergence rate becomes increasingly important as the mesh is 
refined since the log of the spectral radius for single-grid methods generally increases linearly with the mesh 
size, thus making computations on very fine meshes impractical. 

One method that has been successful in accelerating the convergence rate of elliptic problems, i.e., attaining 
a spectral radius independent of the mesh spacing, is the multigrid method (refs. 10 and 11). Although most 
of the existing theory on multigrid methods pertains spe. ific.ally to elliptic equations, it has been shown in a 
number of references (refs. 12 to 22) that the multigrid method can greatly accelerate the convergence rate of 
numerical schemes used for solving the Euler equations. 

One of the earliest applications of multiple grids in sob ing the Euler equations was presented by Ni (ref. 12) 
who used coarse grids to propagate corrections rapidly "throughout the domain. His original idea was first 
incorporated into a one-step Lax-Wendroff method and was later extended for use into predictor-corrector 
types of methods by Johnson (ref. 13). Johnson and ('hima (refs. 14, 15, and 16) subsequently used the 



method to calculate both inviscid and viscous flows over several two-dimensional geometries. In 1984, Mulder 
applied a linear multigrid method to the Euler equations in two space dimensions by using upwind differencing 
to calculate flow over a circular arc and for a weakly barred galaxy (ref. 17). Jespersen also used upwind 
differencing in two spatial dimensions to calculate flow over airfoils (ref. 18). In this approach, the Euler 
equations were solved by Newton iteration where the linear system arising at each step was solved using the 
multigrid method. One of the first uses of the nonlinear multigrid method in accelerating the convergence rate 
for both the two- and three-dimensional Euler equations was reported by Jameson who used central differencing 
in a four-stage Runge-Kutta algorithm to advance the solution (refs. 19 and 20). In two dimensions, recent 
work by Jameson and Yoon also used central differencing and incorporated the multigrid algorithm into some 
implicit schemes with good success (refs. 21 and 22). 

The purpose of the current investigation is to combine the full- approximation-scheme (FAS) multigrid 
method with flux-vector splitting to obtain efficient solutions to the Euler equations in three dimensions. The 
full-approximation scheme for a general nonlinear problem is discussed as well as its implementation for the 
Euler equations. Both V- and W-type cycling strategies are also investigated. Several smoothing algorithms 
are given involving mostly vectorizable computations on the Control Data Corporation VPS-32 supercomputer 
(a CYBER 205 with 32 million words of memory) at the NASA Langley Research Center. In addition, both the 
Steger- Warming and Van Leer splittings are considered for splitting the flux vectors into upwind components. 
Numerical pressure distributions are compared with available experiment for subsonic and transonic flows over 
the ONERA M6 wing and for supersonic flow over an analytically defined forebody. 


Symbols 

A, B, C 

A 

a 

CFL 

c-p 

c 

C/ 

e 

F,G,H 

/energy 

/mass 

G j\i 
grad( ) 

I 


t 

ij 

J 


flux Jacobians in Cartesian coordinates, ^F/5Q, ^G/5Q, and dH/5Q, 
respectively 

flux Jacobians in generalized coordinates, 3F^/^Q, 5G^/(9Q, and 
respectively 

matrix from similarity transformation 
speed of sound 

Courant-Friedrichs-Lewy number 
pressure coefficient 
airfoil chord 
lift coefficient 

total energy per unit volume 

flux vectors of mass, momentum, and energy 

energy component of flux vector 

mass component of flux vector 

grid level N 

gradient operator 

identity matrix 

restriction operator used for transferring functions on grid i to grid z — 1 
prolongation operator used for transferring functions on grid z — 1 to grid z 
restriction operator for the residual 

= x/=I 

X and y components of a unit vector 
transformation Jacobian 
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denotes r/, or ^ 
generic nonlinear operator 
Mach number 

local Mach number in ^-direction 

implicit operator 

unit normal vector 

forcing function added to residual 

pressure 

conserved variables representing m iss, momentum, and total energy per unit 
volume 

change in dependent variables, ^ 
approximate value of conserved variables on grid i 
residual 

Riemann invariant 

components of a vector in x-, and ^-directions, respectively 

forcing function 

entropy 

length of cell face 

arc length, as used in appendix C 

similarity transformation matrix v hose columns are the right eigenvectors of 
the flux Jacobians; also denotes a rotation matrix in appendix B 

similarity transformation matrix whose rows are the left eigenvectors of the 
flux Jacobians; also denotes invert e of rotation matrbc m appendix B 


time 


txi ty , tz 

components of a vector in x-, y-, .ind ^-directions, respectively 

At 

time step 

U,V,W 

contravariant velocities; V" also denotes the volume of a computational cell 

where indicated 

Uo 

a constant vector 

u, u, w 

Cartesian velocities 

Vi 

correction on grid level i 

V 

eigenvector for Fourier analysis 


Cartesian coordinates 

a 

angle of attack 

^ Ps 5 Is 

coefficients for stability analysis 

1 

ratio of specific heats 

A 

incremental change 

1 

difference operator in x-, y-, and ^-directions 


Vs 

A 

A 

Ai, A5 

^xt ^y, Cz 
P 

T 

II 

Subscripts: 

A 

a, b, c, d 

i,j,k 

r 

^,y,z 

Superscripts: 

c 

n 

± 

n 

n 


percent of spanwise location on wing, as used in figures 15, 16, 18, and 19 

diagonal matrix of eigenvalues 

amplification factor 

eigenvalues of Jacobian matrices 

smoothing factor 

general curvilinear coordinates 

components of a unit vector normal to a ^ = Constant face 

density 

time 

relative truncation error 
magnitude of a vector 

large mesh cell in appendix C 

smaller mesh cells in appendix C; b also denotes values on boundary in 
equations (34) to (36) ^ 

cell indices 

discretization on grid N — 
reference values 
spatial derivatives 


most current value of a quantity 
time level 


and eigenvalue contributions; also denotes forward 

and backward spatial differencing or extrapolation where indicated 


quantities in locally orthogonal coordinates 

fluxes and flux Jacobians in generalized coordinates 
symbol where indicated 


; also denotes Fourier 


Euler Solution Method 


Euler Equations in Generalized Coordinates 

t *- dy.«nics, i.e.. the Euler equatione, 

orerr^rua, force,. T.e co„,.r.a.l„„ 


where 


^ ^ ^ 
dr ^ ^ ^ ^ ^ 


( 1 ) 


( P 



1 

J 


pu 
pv i 
pw 


( 2 ) 


4 


e ) 




p{7u 4 ixP 
pUv-\ iyP 
pUw H izP 
I {e + p)U ) 
pV 

pVu + T}xP 
pVv - riyp ^ 
pVw 4 - r)zp 

I {e + p)V J 

pW 

pWu -qxP 


H = — ^ pWv ^qyP 
u 


( 3 ) 




(4) 


(5) 


pW w f (^zP 
. [e->rp)W 

The pressure p is related to the conserved variables through the ideal gas law 

p = (-7 - l)[c - p{v? 4- + w^)/2] (6) 

The equations have been generalized from Cartesian coordinates using a steady transformation of the type 


^ = ri = T]{x,y,z) c = c:{x,y,z) 

T = t 

(7) 

where the contravariant velocity components are 



U = ^j.tt -|- < -\- ^Z'W 


(8a) 

V — r)xU + riy!> + T]zW 


(8b) 

W ^<;xu + ^y ' + <;zw 


(8c) 


The transformation to generalized coordinates is given in appendix A. ^ , r w 

The equations, although written in generalized coordinates, are used m a finite-volume formulation. 
Equation (1) can be interpreted as describing the balance of mass, momentum, and energy over an arbitrary 
control volume. In this connection, the vectors grad(0/^. grad(i 7 )/J, and grad(c)/J represent directed are^ 
of cell interfaces in the contravariant rj-, and <;-direc tions; i.e., m directions along normals to the f 
Constant, rj = Constant, and ? = Constant surfaces, respectively. The Jacobian J represents the inverse of e 
cell volume. Likewise, the quantities pU/J,pV/J, and pW/J represent the contravariant mass flux crossing 
the cell interfaces in the contravariant t]-, and c-directions. 

Flux-Vector Splitting 

The upwind differencing in the present work is effected through the technique of flux-vector splitting. The 
generalized fluxes F,G, and H are split into forward and Ijackward contributions according to the signs of the 
eigenvalues of the Jacobian matrices and are differenced accordingly. For example, the flux m the ^-direction 
can be differenced as 


= 6“F+ f «5^F 


(9) 
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since F has all nonnegative eigenvalues and F has all nonpositive eigenvalues. For the current study, two 
methods are considered for splitting the flux vectors into upwind components. Although the details of each 
method can be found in references 3, 4, 5, 7, and 9, both methods are briefly discussed below. 

The first method presented is the technique outlined by Steger and Warming in reference 3. Since the 
flux vectors are homogeneous functions of degree one in Q, they can be expressed in terms of their Jacobian 
matrices. For example, considering the flux vector in the ^-direction allows F to be written as 


F = AQ = 

dQ 


Using a similarity transformation allows equation (10) to be rewritten as 


(10) 


F = AQ = TAT“^Q 


The matrix A is a diagonal matrix composed of the eigenvalues of A and is given by 


where 


A = 


rAi 

0 

0 

0 

0 


0 

A2 

0 

0 

0 


0 0 0 ■ 

0 0 0 

A3 0 0 

0 A4 0 

0 0 A5. 


Aj 2,3 = u — + ^yV + ' 

A4 = C7 + |grad(^)|a > 

A5 = - |grad(^)|a 

The eigenvalues can then be decomposed into nonnegative and nonpositive components 


( 11 ) 


( 12 ) 


(13) 



II 

1 

(14) 

where 

^ 



* “ 2 

(15) 

Similarly, the eigenvalue matrix A 

can be decomposed into 



A = A+ + A“ 

(16) 


where A+ is made up of the nonnegative contributions of and A" is constructed of the nonpositive 
contributions of A^ . This splitting of the eigenvalue matrix, combined with equation (11), allows the flux 
vector F to be rewritten as 


F = T(A+ + A )T = (A+ + A~)Q = F+ -f-F" 


(17) 


The flux vector F has three distinct eigenvalues given by equation (13) and can therefore be 
sum of three subvectors, each of which has a distinct eigenvalue as a coefficient (ref. 7): 


written as a 
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F = Fi + F2 + F3 


(18) 


where 


P 


pu 

F,=A,l^| p. 

pw 

H + w"^) 

P 

pii i po^^x 

^ pi ± pa^y 

pv ± pa\z 

.«+P±lii^ j 

and the direction cosines of the directed interface in the ^-d rection are 


F2,3 = A4,5J^ 


lx = ^i/|grad(OI 
ly = ^y/|grad(OI 

Iz = ^z/lgradfOl 


( 19 ) 


( 20 ) 


(21a) 

(21b) 

(21c) 


The forward and backward flux vectors F+ and F are formed from equations (18), (19), and (20) by inserting 

A,- = At and Ai = A", respectively. _ 

For supersonic and sonic flow in the ^-direction, i.e., \M^\ = \u/a\ > 1, where ii = t//|grad(OI represents 
the velocity normal to a ^ = Constant face, it should be noi.ed that the fluxes in this direction become 


0 
II 

1 

11 

+ 

(M^ > 1) 

II 

i 

o 

11 

+ 

lA 

1 


(22a) 

(22b) 


The split fluxes in the other two directions are easily obtained by interchanging rj or c in place of 

The fluxes split in the aforementioned manner are not continuously differentiable at zeros of the eigenvalues 
(i.e., sonic and stagnation points). (See ref. 23.) This is illustrated in figure 1 where the split mass flux 
contributions for the one-dimensional Euler equations, nondimensionalized by pa, are shown as a function o 


Mass flux 



Figure 1. Variation of Steger and Warming ??plit mass flux with Mach number. 


7 



Mach number. The gradient discontinuities in the split fluxes are evident as the eigenvalues pass through zero. 
The lack of differentiability of the split fluxes has been shown in some cases to cause small oscillations at sonic 
points that are rarely noticeable for most aerodynamic applications. 

It should also be noted that the Jacobian matrices of F~^ and F that are required for proper linearization 
for an implicit scheme do not have the same eigenvalues as A+ and A“ defined in equation (17). (See ref. 23.) 
However, the Jacobian matrices of F+ and F“ do have the same sign as A+ and A~ so that upwind differencing 
the spatial derivatives remains appropriate (ref. 23). Although A^ are easier to form, their use in implicit 

schemes, instead of the correct linearizations A^, has been shown in many cases to cause severe time-step 
limitations (refs. 3, 24, and 25). 

In 1Q82 a new method of splitting the flux vector was proposed by Van Leer (ref. 5). Here, the fluxes were 
split so that the forward and backward flux contributions blended smoothly at eigenvalue sign changes, i.e. 
near sonic and stagnation points. Just as for the Steger- Warming splitting, it was required that the Jacobian 
matrices 3F“*'/5Q have nonnegative eigenvalues and d¥~ fdC^ have nonpositive eigenvalues so that upwind 
differencing could be used for the spatial derivatives. In addition, it was required that both Jacobians have 
one zero eigenvalue for subsonic Mach numbers which leads to steady transonic shock structures with only 
two transition zones (ref. 5). In practice, when second-order spatial differencing is used, shocks with only one 
interior zone are usually obtained (ref. 9). This feature is not observed with the Steger- Warming flux splitting. 

The three-dimensional splittings of Van Leer were originally given for Cartesian coordinates. The extension 
to generalized coordinates is given in appendix B with the resulting split fluxes given below. Only the splitting 
for the flux in the ^-direction is given, as the others can be obtained similarly. The flux vector F is split according 
to the contravariant Mach number in the ^-direction, defined previously as Afe — u/ct. For supersonic flow 

(|M^I > 1), 

F+ = F F~ == 0 (Aff > -1-1) (23a) 

F“ = F F+ = 0 (Afj < -1) (23b) 

and for subsonic flow (|Af^| < l), 

f /± 

•/mass 

/mass ^ 

^i^jgra^^ /^ass{[ey(-w±2a)/7] + u} 

/mass {[^^(-^ ± 2a)/7] + Mjj 

/± 

./energy 

where 

/mass “ ih 1)^/4 (24b) 

/energy ~ /mass |[~(7 ~ i 2(7 — l)ua + 2a^]/(7^ — 1) -f -j- W^)/2^ (24c) 

For forming F , the direction cosines and \z s-rc given by equations (21) and u is the velocity normal 

to a ^ Constant face. The fluxes in the other two directions are easily formed by interchanging ^ with tj 
or c. In figure 2 the nondimensionalized mass flux using the Van Leer splitting is shown as a function of Mach 
number for the one-dimensional Euler equations. The split fluxes are continuously differentiable at sonic and 
stagnation points; the improvement over the Steger- Warming splitting is apparent. 

Baseline Solution Algorithm 

The baseline algorithm for updating the steady Euler equations stems from a backward Euler time 
integration of the unsteady equations, which yields (ref. 8) 

[l -h At{6- A+ + 6^ A-) + At{6-B+ + 6+B~) + At{S-C+ + 6+C-)j AQ = -At R” (25a) 
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Mass flux 
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Figure 2. Variation of Van Leer split mass flux with Mach number. 


where the residual at time level n is given by 

R" = + 6~G~^ +6^G~ + (5'H+ + (25b) 

The split-flux differences in equations (25) are implemimted as a flux balance across a cell, corresponding 
to MUSCL-type differencing (Monotone Upstream-Centered Schemes for Conservation Laws). (See ref. 26.) 
For example, the flux balance in the ^direction across a cell centered at point {i,j,k) can be written as 

^-F++^+F- = [F+(Q-) + F-(Q+)l.+ i/2)-[F'"(Q")+f'"(Q'^)]r-(l/2) (26) 

The notation F+(Q"’),^.(i/2) denotes the forward flux evaluated using the metric terms at the cell interface 
i + (1/2) and the conserved state variables on the upwind side of the interface, obtained by a fully upwind 
second-order state variable interpolation: 

As seen in figure 3, denotes the average value of Q in the cell centered on Ui,rij,^k) time for 

simplicity, wherever the index notation is {i,j, k) or n, it .s most ofteii dropped. at 2 . -i. 

In equations (25), if second-order differencing is use- 1 on both sides of the equation, Newton iteration 
for the steady Euler equations is obtained as At approaches infinity. The solution however requires a large 
banded block matrix to be solved at each step, a procedure that is generally not feasible because of the large 
number of operations required to invert the system. E/en if the differencing on the left-hand side of the 
equation is reduced to first order, which would not affect the second-order accuracy of the final solution, the 
resulting system of equations usually remains uneconomical to solve. Therefore, the solution is obtained using 
approximate factorization, which splits the implicit operator into a sequence of easily invertible equations. 

When using flux- vector splitting, there are numerous w ays of factoring the implicit operator into a sequence 
of simpler operators (ref. 3). For the results shown below, three ways of factoring are considered. Each 
of the schemes uses first-order spatial differencing on he implicit side of the equation, whereas second- 
order differencing is maintained for the residual calculations. Since the steady state does not depend on 
the differencing of the left-hand side, the final steady-state result will have spatially second-order accuracy 
The computational modules for each of the schemes is shov/n in figure 4. All the schemes employ simple explicit 
boundary conditions. Since only steady-state solutions a e sought, each cell is advanced at its own time step 
corresponding to a given CFL number defined by 

CFL = At {\U\ + |y I + |W1 -t- o [|gr.id(01 + |grad(r?)l -I- lgrad(c)U} (28) 
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O Time level n 
• Time level n + 1 
O Time level n + 1 

(previously calculated) 



(a) Three-factor spatially split. 
( sweep. 



(b) Two-factor eigenvalue split. 
Positive eigenvalue sweep. 


(c) Two-factor combination split. 
e/C+ sweep. 

Figure 4. Computational modules. 

The first scheme considered, is a spatially split algorithm given by 

[I + At(S~A+ + <5+A-)][I + At(S~B^ + + At(S~C+ + S+C~)]AQ = -At R" (29) 

The computational module for the left-hand side of equation (29) is shown in figure 4(a) for the ^ sweep. Since 
he solution at each point is directly coupled to the two neighboring points, the scheme requires the solution 
of a system of block tridiagonals. Similarly, the other two factors also require a block tridiagonal inversion. 

his scheme has the advantage however of being completely vectorizable, and viscous effects can be easily 
inco^orated into the left-hand side. Since the VPS-32 computer is much faster for long vector lengths than 
for short ones the computations in the present implementation take advantage of the large memory available 
on the VPS-32 and solve the block matrix equations over multiple planes simultaneously, thus yielding longer 
vector engths and faster processing rates. As seen in figure 5, the block tridiagonal matrices can be solved with 
vector lengths corresponding to the number of lines in a plane times the number of planes taken. The residual 
calculations, on the other hand, can be made with vector lengths corresponding to the number of points in 
the grid. To decompose the implicit operator into lower (L) and upper (U) matrices (LU decomposition) and 
perform the back-substitutions requires 695 multiplications and additions per factor resulting in a total of 2085 
operations for each sweep through the grid. 

The second method considered to factor the left-hand side of equation (25a) is a two-factor method in which 
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Vectoriza on direction 
for implicit ^ sweep 



Figure 5. Vectorization of LU decomposition sad back-substitution over multiple planes. 

The second method considered to factor the left-hand side of equation (25a) is a two-factor method in which 
the implicit operator is split such that one operator contains the Jacobians with all positive eigenvalues an 
the other operator contains the Jacobians with all negative eigenvalues. The scheme can be written as 


[I-l- -t-(5^B+ -1-6+C )]AQ - At R" 


(30) 


and requires only the solution of block lower triangular equations. Each factor is solved by startup at one 
corner of the grid and solving for each point by marching across the field to the opposite corner. FVom the 
computational module shown in figure 4(b) for the first factor in equation (30) the solution at the center node 
requires that the solution of each of the points behind it be previously obtained. These terms are taken to the 
right-hand side of the equation and added to the residual. A similar procedure is carried out for the second 
factor. For the present implementation of this scheme only 270 operations are required for each factor to 
invert the left-hand side and perform the back-substitul ions, and most of the algorithm is vectonzable. This 
scheme requires that a 5 X 5 matrix be inverted at ea< h point in the grid. These inversions can be carried 
out with vector lengths corresponding to the number of points in the grid. However, since the solution of each 
plane requires that the solution of the previous plane be known, the back-substitution cannot be perforined 
over multiple planes. The maximum vector lengths in the back-substitution process correspond to the number 
of points in a plane. Scalar computations are used to perform the back-substitution along a Ime and comprise 
roughly one-third of the total operations. Although the scheme requires only about 25 percent of the operations 
required for the spatially split scheme, the scalar comput-itions degrade the processing rate significantly so that 
th^ overall computer time per iteration on the VPS-32 is about twice as much as that for the spatially split 

The last scheme considered is another two-factor scheme that is spatially split in two directions, with the 
third direction split according to the sign of its eigemalues. The resulting scheme, which is referred to as 
combination splitting, is given by 


[I -f- At(^^ A”*” 4“ <5^ A 


+ 6-C+)][l + At{ 


•,-B+-h^+B- 


+ 6^C' 


-)]AQ = -At R" 


(31) 
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and the computational module for the first factor is shown in figure 4(c). This scheme also requires the solution 
of block tridiagonal systems and requires about 695 multiplications and divisions for each factor. This scheme 
is completely vectorizable; however, as in the previous two-factor scheme, the solution of each plane requires 
that the solution of the previous plane be known, thereby eliminating the possibility of extending the vector 
operations over several planes. The result is that even though this scheme requires only two-thirds of the 
operations of the three-factor scheme, the computational rate is actually degraded by about 10 percent. A 
summary of operations required to solve the left-hand side for each of the three schemes is given in table I. 

Table I. Operation Counts for Solving the Implicit Operators 


Scheme 

Operations required per point 
per factor with associated 
vector lengths® for — 

Total operations 
per point 

LU decomposition 

Back substitution 

Three-fax: tor spatially split 

550 ML 

145 ML 

695 X 3 = 2085 

Two-factor combination split 

1 

"550 ML 
50 L 


145 L 

745 X 2 = 1490 

Two- factor eigenvalue split 


f75 MP^ 

I5OP 

50 L 1 
,50 S J 


45 s 

270 X 2 = 540 


— vector length corresponds to number of lines in a plane. 

LU— lower and upper. 

ML vector length corresponds to number of lines in a plane times number of planes. 
MP vector length corresponds to number of points in a plane times number of planes. 
P^vector length corresponds to number of points in a plane. 

S — scalar computation. 


Boundary Conditions 

The boundary conditions for the solutions presented below are applied explicitly. On the body, 
the normal velocity is set to zero, whereas the pressure and density are determined by zeroeth-order 
extrapolation from the interior. For subsonic flow in the far field, the velocity normal to the boundary 
and the speed of sound are obtained from two locally one-dimensional Riemann invariants given by 


= u± 


2a 

1 - 1 


(32) 


These invariants are considered constant along characteristics defined normal to the outer boundary 
given by 


/ dx\^ 

[mJ 


(33) 


For subsonic conditions at the boundary, R~ can be evaluated locally from free-stream conditions outside 
the computational domain and is evaluated locally from the interior of the domain. The local normal 
velocity and speed of sound on the boundary are calculated using the Riemann invariants as 


Ub = -{R+ + R-) 
06- ^(R+-R-) 


(34a) 

(34b) 
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The Cartesian velocities are determined on the ( uter boundary by decomposing the normal and 
tangential velocity vectors into components yielding 


Uj, = Ur + kxi'h ~ '“r) 

Vb = + kyi'ib - ^r) ' 

Wb = Wr + kzi Ub — Ur) , 


(35) 


where the subscript r represents values obtained from one point outside the domain for inflow and from 
one point inside the domain for outflow. 

The entropy S is determined using the value from outside or inside the domain, depending on whether 
the boundary is an inflow or outflow boundary. Once the entropy is known, the density on the far held 
boundary is calculated from the entropy and speed of sound as 


Pb = 



(36) 


The energy is then calculated from the equation of st ite. , , . r 

For supersonic free-stream conditions along inflov boundaries, quantities are extrapolated from the 
exterior; along outflow boundaries, quantities are extrapolated from the interior of the computational 

domain. 


Stability Analysis 

In order to examine the stability characteristics 
algorithms considered previously, a Fourier analysis 
in Cartesian coordinates (refs. 8 and 25). Because 
equations and the fact that the three-dimensional 1 
system of convection equations, stability analysis of 
determine stability properties of the three schemes, 
split scheme that would reduce to only a one-factor s 
equations can be written as 

N AQ = -L 


>f the three-dimensional approximate-factorization 
IS conducted on the complete system of equations 
)f the mixed signs of the eigenvalues of the Euler 
juler equations cannot be diagonalized to yield a 
the scalar convection equation is not sufficient to 
(For example, consider the two- factor eigenvalue 
cheme for scalar analysis.) The complete system of 






where, for Cartesian coordinates, 

R" = F+ + 6^F~ + 6yG^ + 6+G“ + (38) 

and N is an implicit operator corresponding to the s heme considered. Linearizing the residual R" as 

R" = Q" + C +67 Q" + C~6+Ct (39) 

and assuming that the Jacobians are locally constan allows the stability to be analyzed by letting 

Q« = A"Uo exp(i^sx] expinsV) exp{iasz) (40) 


where Uq is an initial constant vector. Upon substitution into equation (37) and dividing out the common 
factors, the generalized eigenvalue problem for A, which is the vector of amplification factors, can be 

(N-L)v = NAv (41) 

where N and L are the Fourier symbols of N and L, respectively. The stabilHy characteristics are 
determined by cycling through a fixed number of eacf of the spatial frequencies, in this case 16 frequencies, 

in the range . ^ o 

0 < Ax, 7 Ly, a Az <2 tt 
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for a series of CFL numbers between 0.1 and 50. Each time the generalized eigenvalue problem is 
solved using a routine from the International Mathematics and Statistics Library (IMSL). Each time the 

maximum eigenvalue, the average eigenvalue, and the smoothing factor are determined. The smoothing 
factor is defined as ^ 

n = max(|A|) (42a) 

when 

7t/2 < max(/3s Ai, 7^ Aj/, a* d^z) < 37t/2 (42b) 

This corresponds to the damping of the high frequencies and serves as an indication of how effectively 
the multigrid procedure can accelerate convergence for a given scheme. 

Results are shown using the Van Leer splitting for each of the factorization schemes given previously. 
Identical cases were run using the Steger- Warming splittings with little change in the results, and these are 
therefore not shown. Each result was obtained by using first-order differencing on the implicit side of the 
equation and fully upwind, second-order differencing for the residual computations. All the calculations 
assume Cartesian coordinates, a Mach number of 0.8, and 0° yaw and angle of attack. 

The average eigenvalue, the smoothing factor, and the maximum eigenvalue are shown in figure 6 for 
the three schemes given previously. For the three-factor, spatially split scheme shown in figure 6(a) the 
maximum eigenvalue indicates that this scheme is conditionally stable with a maximum CFL number 



A Maximum eigenvalue 
□ Average eigenvalue 



(a) Three-factor spatially split. 


(b) Two-factor eigenvalue split. 



(c) Two-factor combination split. 


Figure 6. Stability analysis of three-dimensional approximate factorization schemes. Moo = 0.8; a - 0"^ 


of approximately 20. Note that this is contrary to the unconditional 
differencing is used with the same algorithm (ref. 27). The minimum 


instability obtained when central 
smoothing factor occurs at a CFL 
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of about 5 which is somewhat less than the CFL numb ;r where the maximum eigenvalue is lowest. In 
contrast to the spatially split algorithm, both of the tw )-factor schemes shown in figures 6(b) and 6(c) 
appear to be stable for all CFL numbers considered. T ie maximum eigenvalues and smoothing factors 
also exhibit less sensitivity to the CFL number than th ? three-factor scheme with minimum smoothing 
factors also occurring at a CFL number of about 5, However, the two-factor eigenvalue split scheme 
has a somewhat higher smoothing factor than the oth< r two schemes. This indicates that the three- 
factor spatially split and the two-factor combination spl t algorithms are more appropriate for multigrid 
applications than the two-factor eigenvalue split scheme. 

Multigrid Algorithms 

General Algorithm 

The multigrid method used in the current study is the full-approximation scheme (FAS) that appears 
in many references (refs. 10, 28, 29, and 30) and is summarized below. It is most easily understood by 
first considering the solution of a general nonlinear equation 

L(Q) = S (43) 

Equation (43) is to be solved numerically by dividing the domain into discrete cells yielding a system 
of equations to be solved simultaneously at each point a-s 

Ln(Qa) = Siv 

where is the exact solution to the discretized system and Lpf is the discrete analog of the operator 
L. If initial conditions are close enough to the final soiution, equation (44) could be solved iteratively 
by using Newton iteration. This approach however may be prohibitively expensive if the number of 
unknowns is large as typically occurs in multidimension.il problems. Many other iterative schemes have 
therefore been devised that require significantly fewer operations. After a few iterations, however, these 
methods generally exhibit a slow convergence rate, thui; reducing the residuals by a very small amount 
each time (ref. 29). The reason for the slow asymptotu convergence rate is inadequate damping of the 
low-frequency errors (ref. 11). 

The multigrid method efficiently damps the low-1 requency errors by using a sequence of grids 
Go,Gi, ...,Gpi, where G^ denotes the finest grid from v.hich successively coarser grids can be formed by 
deleting every other mesh line in each direction. In this context, the high-frequency-error components 
on a given grid are those that cannot be resolved on th« next coarser mesh because of the increased grid 
spacing. If an iterative method is chosen that quickly clamps the high-frequency errors on a given grid, 
then the remaining errors will be only the low-frequency smooth components after a few iterations. A 
sequence of coarser grids can then be used to accelerate the convergence rate on the finest grid by reducing 
the remaining low-frequency errors since some of these i ame frequencies appear as high-frequency errors 
on a coarser grid. Therefore, the errors on the fine grid that are usually responsible for slow convergence 
are quickly damped when using the coarser grids where the computations are relatively inexpensive. 

In order to use the coarser grids, it is necessary to obtain an equation on the fine mesh that can be 
accurately represented by the coarser mesh. It is important to realize first that neither the solution nor 
the high-frequency-error components on the fine grid c.in generally be resolved on a coarser grid. The 
high-frequency errors however can be sufficiently damped on a fine grid by using appropriate iterative 
schemes so that the remaining errors will be composed of only smooth, low-frequency components that 
can be adequately represented on coarser meshes. For his reason, it is necessary to obtain an equation 
on the fine mesh in terms of the errors. 

When solving in an iterative fashion, equation (44) is solved approximately at each step as 

~ ' 4 " ( 4 ^) 

where is the most current approximation to Q;v and is the residual that will be zero only when 
= Qpi and, hence, the exact solution is obtained. Subtracting equation (45) from equation (44) yields 
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an equation on the finest grid in terms of the residual 


^n{Qn) ~ ^n{^%) — ~R-N (46) 

If the high-frequency errors have been previously smoothed, then the fine-grid residual equation (46) can 
be adequately approximated on a coarser mesh by 

L;v-i(Q7V-i) = (47) 

where ^ and ^ are restriction operators for transferring the dependent variables and the residual 

from the fine grid to the coarse grid, respectively. Here, serves as an initial approximation to 

the solution on the coarse mesh, whereas Q^v-i is the exact solution which is the sum of the initial 
approximation and a correction (ref. 30). Since the full solution is computed and stored on each grid level 
as opposed to only the corrections, this process is referred to as the full- approximation scheme (FAS). 

On a sufficiently coarse grid, equation (47) can be solved exactly using a variety of numerical techniques 
to obtain Qtv_i from which the correction can be formed as 


Vn-i = Qat-1 - 1% (48) 

This can then be passed up to the fine grid and used as a correction to q^, which is replaced by its 
previous value plus the correction 

ijv-1 Vat-1 (49) 

This process yields a simple FAS two-level algorithm where the operations on the coarse grid (eqs. (47) to 
(49)) that are used to update the fine-grid solution are termed the coarse-grid correction. Often, however, 
the exact solution of equation (47) can still be quite expensive to obtain. Also, since the correction 
on the coarse grid serves only as an approximation to the fine-grid correction, the exact solution of 
equation (47) is not required. Therefore, instead of solving equation (47) to completion, several iterations 
can be carried out to get a reasonable approximation to Qat_i. After each iteration of equation (47), the 
equation satisfied by q/^_, is given by 


LA'^-l(q^_l) — ^(— Rat) + LAT-i(lj^ ^q^) + R-AT-l (50) 

which differs from the solution of equation (47) only by the residual term Rat-i which will be zero when 
q^_l = Qat- 1- If the errors are smooth, then subtraction of equation (50) from equation (47) yields 
an equation that can be well-represented on yet a coarser mesh, Gat- 2- Writing this equation on Gat _2 
yields 

^N-2{Qn~2) = h'AT-2(ljv-l*l^-l) fjv-i(~H-7V-l) (51) 

where equation (50) is used to determine Rat-i- The solution may be obtained by solving equation (51) 
exactly, by approximating by several iterations, or by introducing more coarse-grid levels. On all coarse 
grids, one or more FAS cycles (smoothing followed by coarse-grid correction) are completed. In this 
manner, each of the coarse meshes is used to obtain a correction for the solution on the next finer mesh. 
Since only the equations for smooth error components may be represented well on coarser grids, it is 
important to only pass corrections and not the full solution from a coarse grid up to the next finer one 
(ref. 28). 

Using equation (45), note that equation (47) can be recast 

LAT-i(QAr-i) = Sat_i - brAT_i == Pat-1 (52) 

where 


16 


s^v-i = 

tat-i = 


(53) 

(54) 


Here, r;v_i is the relative truncation error (or defect correction) between the grids so that the solution on 
the coarse grid is driven by the fine grid, and the defect correction accounts for the difference in truncation 
error between the coarse and fine grids (ref. 28). The a ialogous equation for grid G;v-2 is given by 

L7V-2(QiV-2) = S.v-2 + ^N-2 


where 


c - 
^N-2 — 




(56) 


^N-2 = “ ijv-; 

Note that the relative truncation error on the N — 2 grid is the sum of the relative truncation error 
between grids N and N — 1, as well as N — 1 and N — 2. 

Algorithm for Euler Equations 

For the steady Euler equations in generalized coordinates, equation (44) can be written as 


Liv(Qiv) = H 6+G- + 5-H+ + = 0 (58) 

In the multigrid-solution process, a forcing function arises on the coarse grids from restricting the residual 
equation on a fine mesh down to the coarser one. Since S = 0 in equation (44) for the Euler equations, 
the forcing function is the relative truncation error between the two meshes. The resulting equation to 
be solved on any mesh G, can be written as 

L,(Q.) = r, (59) 


where Tj = 0 on the finest mesh and is the relative truncation error on each of the coarser grids. The 
solution of equation (59) is updated by introducing a l ime derivative of the dependent variables to the 
left-hand side so that the solution can be advanced in f ime using the approximate-factorization methods 
previously described. The resulting scheme written on mesh is given by 

N Aq? = - At [Li(q^; - r^] = -At (60) 

where N is the implicit operator of the scheme considered and Lt (qJ) on the right-hand side is due to the 
linearization of L,(Q,) from the backward Euler time integration. Note that even on the coarse meshes 
where is nonzero, equation (60) maintains the same form as the equation on the fine mesh. The result 
of this is that the coarse meshes can be updated usini' the same scheme as that used on the fine mesh 
with only a slight modification to the right-hand side. 

Several strategies exist for deciding when to switch from one grid level to another, and they generally 
fall under the categories of fixed- or adaptive-cycling algorithms. The strategy used in the present study 
is a fixed-cycling strategy in which each global cycle consists of a set number of FAS cycles on each of the 
coarser grids. Recall that one FAS cycle on any grid consists of a smoothing step followed by a coarse-gnd 
correction. A predetermined number of iterations are j.erformed on each grid level to smooth the errors. 

The conserved variables are transferred to the next coarser grids each time by the rule 


Q.-i=ir^Qi 

where is a volume- weighted restriction operator tliat transfers values on the fine grid to the coarser 
one and is defined by 

and the summations are taken over all the fine-grid ceils that make up the coarse-grid cell. As shown in 
appendix C, restriction of the dependent variables in this manner conserves mass, momentum, and energy 
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in each of the cell volumes. The relative truncation error is calculated on the coarse grid as 

r,_i = L,_i(i:-V)-ir'R. (63) 

where 1'“^ is the restriction operator for the residual defined as 


~ (64) 

where, again, the summation is over the cells on the fine grid that make up the coarse-grid cell. By 
summing the residuals, the surface integral of the fluxes crossing the cell boundaries on the coarse grid 
IS the same as would arise by integrating around all the fine-grid cells making up the coarse grid. (See 
appendix C.) Several iterations of the approximate factorization scheme can be conducted to get an 
approximation to the steady solution on with the right-hand side modified to include the relative 

truncation error. If only one coarse grid is used to correct the finest grid, the result is the simple FAS two- 
level cycle. On the other hand, if more grid levels are introduced so that one or more FAS cycles can be 
recursively carried out on each subsequent coarse-grid level to get a better approximation to Qjv_i , then 
a multilevel algorithm results. When only one FAS cycle is carried out for each of the coarser grid’s, the 
resulting global cycling strategy is termed a V-cycle and is depicted in figure 7. Another cycling strategy 
of interest, which is shown in figure 8, is termed a W-cycle and results when two FAS cycles are used on 
each of the coarser meshes. Results will be shown in the next section using both types of cycles where, 
on the coarsest mesh, three smoothing iterations are performed in lieu of solving the equations exactly 
on the coarsest mesh. The corrections on coarse meshes are passed to the next finer mesh using trilinear 
interpolation with no additional iteration steps between meshes. When a W-cycle is used, however, note 
that an iteration is carried out at the beginning of each FAS-cycle correction in order to smooth the high 
frequencies. 
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In order to clarify the multigrid procedure further 
exemplary case where three grid levels are used in a 

1. Start on the finest grid and smooth the errors 

2. Calculate the residual on the fine grid from equ 
side of equation (58) and S;y = 

3. Restrict the dependent variables to the first co 

4. Restrict the residual from the finest grid to G 
truncation error using equation (63). 

5. Calculate the right-hand side of equation (60) a 
to smooth the errors on this grid so that a coarser gi 

6. Calculate the residual on this mesh using equa 


the overall process is summarized as follows for an 
^/-cycle: 

by doing one iteration of equation (60) with = 0. 
ition (45) where LTy(q^) is given on the right-hand 

irse grid Gat-i by using equation (61). 

using equation (64) and calculate the relative 

ad update the solution on mesh Gjv_i. (This serves 

id can be introduced.) 

cion (50). Note that this can be written as 


RjV-l = 


Since r^-i has been previously calculated, the residual is easily calculated by simply calculating 
LjV-l(q/v-i) from the most current values of the dependent variables on the mesh and subtracting 

^ 7. Restrict the dependent variables on Gyv-l to > by using equation (61). 

8. Restrict the residual from equation (65) to the AT - 2 grid and calculate r^v-2 from equation (63). 

9. Calculate the right-hand side of equation (60) and update the solution on this inesh. Since this is 
the coarsest mesh used in the present example, three iterations of equation (60) are undertaken to get an 
approximation to Q^_ 2 - During each step, the right -hand side is updated to use the most current values 
of the dependent variables in Ljv-2(q^-_2)' Note that rjy-2 will not change. 

10. Calculate the correction on this mesh to give 

a/ - ..iC (66) 

Vjv-2 — “ -fjv-iqiV-i ' 

n. Pass the correction to the next finest mesh using trilinear interpolation and update the solution 

Note that steps 5 to 11 make up one FAS cycle on grid N - 1 where steps 6 to 11 constitute a coarse- grid 
correction. At this point, if a W-cycle were being employed, another FAS cycle (steps 5 to 11) would be 

repeated to update q/y_i further. 

12. Calculate the correction on the N — I mesh is 

xr C 

^N-l — *lN -l 

13. Pass this correction to the finest mesh and u pdate the solution to give 

q^- q^ 

14. Perform one smoothing iteration using equa ion (60) to smooth the errors. 

Presentation of Results 


ONERA M6 Wing 

Three-dimensional subsonic- and transonic-flow 
presented. Comparisons are made with experiment 
corresponding to conditions for which viscous effects 
airfoil sections with a planform swept 30° along the 1 
0.56. Solutions are obtained for two mesh types, C- 
around the airfoil profile. The C-H mesh has uni 
C-O mesh wraps around the wingtip, consequent! 


computations over the ONERA M6 wing are now 
il data at a Reynolds number of 11.7 X 10® (ref. 31), 
are relatively small. The wing consists of symmetrical 
i!ading edge, an aspect ratio of 3.8, and a taper ratio of 
H and C-O, both of which are C-type mesh topologies 
form spacing in the spanwise direction, whereas the 
y leading to a more precise definition of the actual 
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rounded-tip geometry tested in the experiment. The C-0 mesh has been generated with a transfinite 
interpolation procedure given in reference 32. The C-H mesh was obtained by simply stacking a series of 
two-dimensional cross sections along the span. The surface meshes for both are shown in figure 9. 

The first computation is on the ONERA M6 wing at transonic conditions; a Mach number of 0.84 and 
an angle of attack of 3.06°. Figure 10 shows the effect of using multigrid on the residual and lift coefficient 
histones for a mesh with over 210000 points. The mesh is a 193 x 33 x 33 C-H mesh. This corresponds 
to 193 points along the airfoil and wake (110 of which are on the airfoil), 33 points approximately normal 
to the airfoil, and 33 points in the spanwise direction (17 of which are on the wing planform). For this 
case the Van Leer splittings are used with a V-cycle and four grid levels (a fine grid and three coarser 
ones). The multigrid method is very effective in accelerating convergence of both the residual and lift 
coefficients. The residual is reduced to “machine zero” in 400 cycles, whereas the single-grid method has 
reduced the residual only between one and two orders of magnitude. The benefit of multigrid is especially 
pronounced in the lift coefficient history where the lift coefiicient value is obtained to within 0.1 percent 
of its final value in 41 cycles. This is a dramatic improvement over the single-grid result that required 
more than 400 iterations to converge to the same level of accuracy for the lift coefficient. It should be 
noted that for all the cases considered, several multigrid V-cycles (usually five) were run with first-order 
spatial differencing before switching to second order. 



(a) 97 X 17 X 17 C-H mesh. 



(b) 97 X 17 X 17 C-O mesh. 

Figure 9. Surface mesh for ONERA M6 wing. 
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Figure 10. Effect of multigrid on 
convergence on ONERA M6 wing 
with 193 X 33 X 33 C-H mesh. 
Moo = 0.84; a = 3.06''. 
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A comparison of convergence rates between the three schemes discussed earlier is shown in figure 11 for 
identical conditions as given above with the exception that only every other point from the 193 x 33 x 33 
mesh is used, thus resulting in a 97 x 17 X 17 C-H mesh. For this mesh size, only two coarser grids are 
used. The three-factor, spatially split algorithm demonstrates a faster rate of convergence than either 
of the two-factor schemes. The spectral radius indicating the error reduction per cycle is 0.898. The 
two-factor scheme in which the implicit operator is split according to the sign of the eigenvalues displays 
the slowest convergence rate with a spectral radius of 0.93. It should be pointed out however that even 
though the spectral radius using this scheme is not as good as that for the spatially split scheme, this 
still represents a good improvement over a corresponding single-grid spectral radius of 0.98. All the runs 
on the 97 X 17 X 17 meshes were made at a CFL number of 7. This was determined experimentally to 
be nearly optimum and agrees well with the CFL number for best smoothing predicted by the stability 
analysis. 

In 64-bit precision on the VPS-32, the computational rate using a three-grid V-cycle for the three-factor 
scheme is about 75 psec per grid per cycle, whereas th( two-factor eigenvalue-split and combination-split 
schemes exhibit computational rates of 140 and 85 /.sec per grid point per cycle, respectively. The 
computational rates obtained using one grid level are decreased about 15 percent from those using three 
grid levels. Computational times are decreased approximately 40 percent when the computations are 
done in 32-bit precision on the VPS-32. Because of th.; higher performance of the three-factor spatially 
split algorithm in both convergence rate and computational rate, it is used exclusively in the results that 

follow. 

A comparison of the residual convergence histories for a W-cycle and the previously used V-cycle is 
shown in figure 12 for the 97 x 17 x 17 C-H mesh. \n improvement in the convergence rate using a 
W-cycle is apparent. In addition, the lift coefficient is obtained to within 0.3 percent of its final value in 
only 14 cycles and to under 0.1 percent in 24 cycles. 1 his is an improvement over the V-cycle that took 
37 cycles to get the error in lift below 0.1 percent. Although the work involved for a W-cycle is more than 
that for the V-cycle because of the extra smoothing iterations on the coarser grids, the time required per 
cycle increased only about 13 percent over a V-cycle. Therefore, even though more work is involved for 
each cycle, a net gain is achieved by employing the W cycle. A summary of results for this case is gi-ron 
in table II for 193 X 33 x 33 and 97 x 17 x 17 grid poinrs for both C-H and C-O types of grids. This table 
includes the spectral radius based on cycles, the numb. T of cycles required to obtain the lift coefficient to 
within 0.3 and 0.1 percent of its final value, and the number of cycles required to obtain the lift to hve 
significant digits. Note that the number of cycles required for the W-cycle to obtain the lift coefficient is 
relatively insensitive to the number of grid points, thus indicating that the multigrid algorithm is working 
successfully. 

Figure 13 shows the upper-surface pressure distributions on the 193 x 33 x 33 C-H mesh as well as on 
the 193 X 33 X 33 C-O mesh. The wing under these conditions exhibits both a swept shock emanating 
from the root leading edge and a nearly normal shock emanating from the root. Both shocks coalesce at 
about 80 percent of the span to form a single shock. F igure 14 shows upper-surface pressure contours for 

the C-O mesh. 
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Two-factor eigenvalue split 

Two-factor combination split 

Three-factor spatially split 



Cycles 


Figure 11, Comparison of convergence rate for 
the three schemes on ONERA M6 wing with 
97 X 17 X 17 C-H mesh. Moo = 0.84; 
a = 3.06°. 



Cycles 

Figure 12. Comparison of convergence rate for 
V-cycle and W-cycle on ONERA M6 wing 
with 97 X 17 X 17 C-H mesh. Moo ~ 0,84; 
a ^ 3.06°. 



(a) 193 X 33 X 33 C-O (b) 193 x 33 x 33 C-H 

mesh. mesh. 


Figure 13. Upper-surface variation of pressure 
coefficient on ONERA M6 wing. Moo = 0.84; 
a = 3.06°. 



Figure 14. Pressure coefficient contours on upper 
surfa.ce on ONERA M6 wing. Moo = 0.84; 
a = 3.06°. 
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Table II. Summary of Results for ONERA Me Wing With Moo 


= 0.84 and a = 3.06^ 



Number of cycles requireci to obtain q within— 

Spectral 

raxlius 

Type of 
cycle 

0.3 percent of 
final value 

0.1 percent of 
final va iue 

Five decimal 

places 1 


97 X 17 X 17 C-II mesh size 



V-cycle 

W-cycle 

20 

14 

37 

24 

75 

42 

0.898 

0.871 


97 X 17 X 17 C-O mesh size 



V-cycle 

W-cycle 

34 

15 

45 

27 

91 

44 

0.912 

0.879 


193 X 33 X 33 C- d mesh size 



V-cycle 

W-cycle 

37 

12 

41 

23 

153 

47 

0.948 

0.923 


193 X 33 X 33 C- 0 mesh size 


V-cycle 

W-cycle 

27 

14 

68 

19 

149 

47 

0.952 

0.926 1 


the agreement with experiment at the m vr^primental data in figure 16. The computations 

fglx 33^31 ?H mesh Tnd C-oZ^h. Figure 17 shows the residual history for both the 97 x 17 x 1 7 C-O 
iy«j X do X 00 XI iiieoii ,. . 1 xi rr^‘^A \F T'Vip ronvercence rate on tne L/'-rl 

iLt Sr 2 t.r 6 .rer ia -a. th«e showa . table .. is ,ve„ In 

table III. 

,. .. ^. fi,.a 07 V 1 7 V 1 7 r-f t mesh and the 193 X 33 x 33 C-H mesh are compared 

The pressure distributions on the 97 L^^j. splittings and in figure 19 using the 

with experiment at six spanwise stations in g g hoth meshes are essentially identical 

Steger-Wamting .plittinp. At a nres contpuL on the 

Steger- Warming and Van Leer splittings give nearly identical results. 
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Figure 17. Multigrid convergence for C-H and C-0 meshes or ONERA M6 wing for V-cycle, Moo = 0.699; q = 3.06°. 
Table III. Summary of Results for ONERA Mfi Wing With Moo = 0.699 and a = 3.06° 

i to obtain ci within — 

ent of Five decimal Spectral 

alue places radius 

I mesh size 




Analytic Forebody 

The last test case considered is an analytically defined forebody for which experimental data are 
available at supersonic Mach numbers (ref. 33), The grid used was a 49 x 49 x 49 grid with a symmetry 
plane along the centerline. Figure 20(a) shows a cross s action of the grid around the body, and figure 20(b) 
illustrates the body surface together with contours of tl le static density in the upper symmetry plane. The 
conditions correspond to a free-stream Mach number of 1.7 and an angle of attack of 0®; this leads to an 
oblique shock at the nose and supersonic flow over the entire length of the body. In this case, it is known 
that a pure axial marching scheme can obtain the solut on quite efficiently, as the equations are hyperbolic 
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Figure 18. Comparison of experiment and inviscid calcula- Figure 19. Comparison of experiment and inviscid calcula- 
tions using Van Leer splittings on ONERA M6 wing for tions using Steger- Warming splittings on ONERA M6 

Moo = 0.699 and a = 3.06^. wing for Moo = 0,699 and a = 3.06°. 




llgStL^n this fully ^supersonic flow. The occurrence of the subsonic pockets presents no difficulty to 
the present approach. 



(a) Cross section of grid around body. 



(b) Grid on body surface with con- 
tours of static density in upper 
symmetry plane. 


Figure 20. Grid for analytic forebody. 


The residual and lift histories obtained using a V cycle and «'>' >“VCcSs‘'rthf !S 
figJe 21 As can be seen, the residual is “ 3 “^ “hie’^rf bl^d on “e 

experimental data over both the lower and upper surfaces of the body. 
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Figure 21. Convergence history for an? lytic forebody. Moo - IT; a 0 
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Figure 22. Comparison of calculated 
Moo = 1.7; a = 0°. 

Concluding Remarks 

A multigrid algorithm has been developed, analyzed, and applied to the three-dimensional flux-split 
Euler equations m generalized body-fitted, curvilinear coordinates. Three implicit smoothing operators 
ere studied, as well as two flux-splitting schemes. Applications were demonstrated for subsonic and 

“re 

A linearized stability analysis of the three-dimensional system of difference equations representing the 

rree..tln, co„diC. S„tee„ 

tbrp f t m each spatial direction were examined. The results indicated that the 

diff!r/n conditionally stable for the upwind, flux-vector-splitting 

difference scheme. This is m contrast to the unconditional instability of the three-factor scher^e with 

un^ond ? eigenvalue-split and combination-split schemes were founTto be 

^^oo^h t ^ However, since extremely large time steps are not necessary for a multigrid 

smwthmg algorithm, the conditional stability of the three-factor scheme is not a penalty ^ 

spaSllv sobri^Z * f “'"“‘Srid smoothing operators, the three-factor, 

spatially split operator was found to be superior to the two-factor operators. Its superiority is due 

to a combination of a better smoothing rate and more complete vectorization 

V-cyde. Either algorithm was an order of magnitude faster than the single-grid iteration. The multigrid 
me^od proved efficient for subsonic, transonic, and fully supersonic flows ^ 

tested, the Steger- Warming and Van Leer splittings, both yielded 
^milar results for pressure distributions. Although not shown, the convergence rates using the Steger- 
Warmmg splittings were very similar to those obtained with the Van Leer splittings ^ ® 


NASA Langley Research Center 
Hampton, Virginia 23665-5225 
August 25, 1988 
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Appendix A 


Transformation to Generalized Coordinates 

The three-dimensional Euler equations in Cartesian coordinates and strong conservation- law form are given 

<9Q ,dG dH 


— 1 I 

dt dx dy dz 


where 


Q = { pu 


pu^ -I- p 


(e + p)u 


G = < pv^ + p 


(e -t- p)v 


pw^ p 
(e -t- p)w 


P= il - ^) 


When using the chain rule and the body-fitted coordinate system given by the steady transformation 

^ = ^{x, y,z) ry = r]{x, y,z) c = <r(i, y,z) r = t (A3) 

the Euler equations can be recast as 

dQ dF dF dF ^ dG dC- dG dU dH dFl , . 

Now, again using the chain rule, the derivatives with respect to r, rj, and ^ can be written m matrix form 
^ / a \ f 9 ^ 

^ r 1 0 I ) 0 “ ^ 

^ ^ 0 Xj if Zf < ^ y (A5) 

^ 0 X^ Ir, Zr, ^ 

d Lo x^ Vi z<;\ a 

I g? J '' 

from which Cramer’s rule can be used to solve for the x, y, z, and t derivatives. Using these to evaluate the 
metric terms gives 

ix = J{Vr)Z<; - Zr)V<i) Vx ~ ~ ~ j 

iy^J{Zr^X^-XrjZ<;) T)y = J {x^Z^ ~ Z^X<;) ^y ^ J{z^Xr, - X^Zrj) > (A6) 

iz = J{xr)y<;-yri^‘i) 'Hz = J{xt;y^ - y,iX^) Cz = j 


where 


wuerc V f A 7\ 

J~^ = X^{yr^Zc^ - Zr^y^;) - y^{xr,z - Z,fX(;) + Z^{Xrfy(,~ yrt^a) V ) 

In order to regain the strong conservation-law form, equa ion (A4) can now be multiplied through by J and 
rearranged using the chain rule on certain terms. For example, the term J can be rewritten as 
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After rewriting all the appropriate terms and noting that many parts of the resulting equation can be shown 
to be zero by substituting equations (A6) for the metrics, the Euler equations can be written in generalized 
coordinates and maintain the strong conservation- law form as 


dQ dF dG dH 

— -I 1 1 

dr dr) 


(A9) 


where 


F = j-i(exF + ^j,G + e^H) 
G = J-^ir]:,F+r)yG + rj^Fl) 
Fi = J ^(?iF -1- CyG CzH) ^ 


(AlO) 


Substituting equations (A2) into 


(AlO) allows the flux vectors to be further written in an alternate form as 



pu ] 


pV 1 


pW 'j 


frnU + (xP 

1 

puV + rjxP 

1 

puW -1- CiP 


pvU + iyp 
PwU + Gp 
. [e + p)U ; 

, G = -. 

pvV + T)yP 
pwV -1- r)zP 
. {e + p)V . 

■ « = 7 ' 

PVW -1- tJyP > 
fnvW -1- CzP 
. {e+p)W . 


where C/, V, and W are the contravariant velocities defined as 


(All) 


U = ^xU -(- ^yV + 

V — r)xU -f- r)yV -I- r)zW * 
W = CxU -I- qyV 


(A12) 
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Appendix B 


Splitting Flux Vectors in Generalized Coordinates 

The Van Leer method of splitting the flux vectors was c riginally given only for a Cartesian coordinate system 
(ref. 5). For example, the flux-split vectors in the i-direc ion were given in terms of the local one-dimensional 
Mach number Mx = u/a. For supersonic flow, i.e., \Mx\ > 1, we have 


F+ = F F" = 0 {Mx> 1) 
F+=0 F"=F (Mi < -1) 


(Bl) 


F± = < 


/mass 


and for subsonic flow, i.e., \Mx\ < 1) have 

±pa[^(Mi:tl)]^ =/= 

/mass [(7 - l)w ± 2 a] 

/ri ass^ 
fnass^ 

/^ass [lil - 1)« ± 2a]2 j [ 2(72 - 1)] -h (u2 + u;2)/2| 


(B2) 


For many applications, however, it is advantageous to construct generalized (body-fitted) coordinate systems 
of the type 

^=i{x,y,z) ri = T]{x,y,z] <; = q[x,y,z) r = t (B3) 

where, in the present work, the transformation is chosen s > that the grid spacing in the computational domain is 
uniform and of unit length. In the discussion that follows, the superscript (*) indicates variables m generalized 
coordinates and the superscript (") indicates variables in a localized Cartesian coordinate system. If no 
superscript is used, Cartesian coordinates are assumed. The strong conservation form of the Euler equations 
in generalized coordinates is given by 


+ + + (B4) 

dr d) d<; 


For the purpose of determining a generalized splittin ; for F, only the derivatives in the ^ and t-directions 
are considered, whereas the rj and <; derivatives are treated as source terms. For determining the splitting 
of F, equation (B4) is transformed by a local rotation matrix in order to decompose the flux vector F into 
components normal and tangential to a ^ = Constant cell face. The rotation matrix is given by 


1 0 0 0 0 ') 

0 ^ 

0 tx ty fl 

0 fx fy Tz 0 

0 0 0 0 1 


(B5) 


where K L. and \z are components of a unit vector h normal to a ^ = Constant line. The f, and r, are 
components of vectors that are normal to n and to each other so that the three vectors form a loca ized 
Cartesian coordinate system. Note that an infinite number of vectors normal to n exist that form a localized 
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Cartesian coordinate system. These vectors however are arbitrary and their exact specification is unnecessary. 
Multiplication of equation (B4) with the matrix T then yields 


Qt + = -TG,, + TtQ + T^F - TH<^ 


(B6) 


where 


Q = tq = - ^ 


F = TF = 


|grad(^)| 


r p 

pu 
pv 
pw 
e 

' pu 
puu + p 
puv 
puw 

I (e + p)u } 


(B7) 


(B8) 


The rotated velocity component u is the velocity normal to a line of constant representing the scaled 
contravariant velocity component, whereas v and w are normal to u and to each other. Thus, 


U = ixU + ^yV + 

V — txU + tyV + izW 
W = fxU + fyV + fzW 


m 

(BIO) 

(Bll) 


The transformed flux F is of the same functional form as the Cartesian flux vector and thus can be split 
according to any splitting developed for Cartesian coordinates. Therefore, equations for both the Steger- 
Warming and Van Leer splittings can be used to split the flux vector F after replacing the Cartesian velocity 
components it, v, and w by the rotated velocity components u, v, and w, respectively. Applying the rotation T 
to equation (B4) simply allows the flux vector to be split in a one-dimensional fashion along a coordinate axis 
perpendicular to the cell interface. After splitting F, the appropriate splitting for F is determined by applying 
the inverse transformation matrix T“^ to equation (B6), thus leading to 


Qt + (F+ + F-)^ -I- G,, -1- = 0 


(B12) 


with 


T~^F^ _ 


|grad(0| 


/mass i 2a)/'yj -f- uj 

fmasa | ^ 

/mass I ± 2a)/'yj ■+■ 

fe 


where 


/energy 


/mass ■“ ^po>{M^ ± j 4 

/energy ^ /mass | “(T ~ i 2(7 — \)ua 2 a^j j ( 7 ^ — 1) j> (B13) 
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Note that the inverse transformation restores the original form of the equations; i.e., no additional source terms 
arise and the form of G and H is unaffected. This allows splitting of G and H similar to the splitting of F 

"^^Dup^cTting the derivation with the Steger- Warming Cartesian splitting, starting with equation (B 8 ), yields 


F = Fi +F 2 + F 3 


(B14) 


where 


- Igrad^- 2^1 
J 7 


lgrad(01 -^4,5 

^2,3 = j ^ 


Ai 


P 

pu 

pv 

pw 

+ W^) 

p 

pu ± pa 
pv 
pw 

€ -f p i pat) 

U 


jgrad(01 |grad(^)| 


where 


A 4,5 ( /±g|grad( 0 | 

“ |grad(01 ” lgrad(OI 


Applying the inverse transformation to equation (B14) gn^es 

F = T'^F = Fi + F 2 + F 3 


7 — 1 


F2,3 = H5J2:; 


P 

pu 
pv 
pw 

I ^ + 1;^ + w"^) } 
p 

pu ± pa^x 
pv i pa^y 
pw ± pa^z 

J 


This result is identical to the generalized splitting given n references 3 and 7. 


(B15) 

(B16) 

(B17) 

(B18) 

(B19) 

(B20) 

(B21) 
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Appendix C 


Restriction Operators 

transferred from a fine mesh to a coarse mesh so that the mass, momentum, 
the r the coarse-grid cell are the same as those contained in the fine-grid cells that compose 

he coarser grid. For a two-dimensional example, referring to figure 23 shows that the mass, momentum and 
energy m any given cell are given by QV, where V is the volume of the cell and Q is the column vector of 
dependent variables representing the conserved quantities per unit volume. Since cell ^ is comprised of the 

tran T \ ^ relationship that conserves mass, momentum, and energy is easily established for 
transferring the dependent variables. Thus, 


o _ + Q&V6 + QcVc -I- 


(Cl) 


writ"' form Tie steady Euler ec„a,i„„s can be 


/ 


F - n ds = 0 


(C2) 


integral is the surface integral over the volume considered. F is the flux of mass, momentum, 
energy across the boundaries, and n is the outward pointing unit vector normal to the boundary. 


or 


2A 


3A 



2c 

1 

1 

1 

2d 


3c 

c 

1 

1c| 3d 

1 

d 

Id 


4c 

1 

t 

1 

4d 



2a 

A 

1 

r 

1 

2b 


3a 

a 

1 

la [ 3b 

1 

b 

1b 


4a 

1 

1 

1 

4b 


— ^ 


4A 




1A 


Figure 23. Control volumes for restriction of dependent variables and residual. 

Considering the two-dimensional case given above for simplicity shows that the integral around the large 
volume shown m figure 23 is given by ^ 


J F nds-{F ■ n)i^SiA + ( F ■ A)2aS2a + ( • 11)3^53^ -|- {W ■ 11)4^54^ 


(C3) 

( F" • n)i4'5i>i = ( F ■ n)i55i{, -f ( F • n)jj5irf 

Also, note that since the outward pointing normals on adjacent cell boundaries point in opposite directions 
several terms that share a common boundary will cancel. For instance 


(F' • n)ic-S'ic — — ( F ■ n)^jSs(i 


(C5) 
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Therefore, by performing the integrations around each of tiie smaller cells and adding them together, it is seen 
that the integral around the larger cell is simply the sum ( f the integrals around each of the smaller ones: 

f F • h ds = [ F • n ds (C6) 

Jl.c. 's.c. 


where the limits L.C. and S.C. indicate large cells and small cells, respectively. 

In order to better relate this specifically to the Euler equations, consider the 
given by 

j pV h d. = 0 


steady continuity equation 


(C7) 


U — ui vj 


(C8) 


r . r '■ grad(A:) 


(C9) 


A normal to a ^ = Constant line is given by using A: = ^ in equation (C9), and the normal to an rj — Constant 
line is obtained by using r) in the same manner. Now, the length of each face is given by 


S,= 


lgrad(fc)| 


(CIO) 


where k is again chosen to be ^ or 77 , depending on which hice is desired, and J represents the inverse of the cell 
volume (i.e., the Jacobian). Using equations (C3) (with F = pU), (C7), (C8), and (C9) allows the integral 
around cell a in figure 23 to be written as 



(p^\ _ (e^\ + 

V -7 / la V *7 /' 3a V <7 / 2a \ -7 / 4a 


(Cll) 


where 

U = 

V = r]xU + qyV 

When calculated numerically, equation (Cll) is simply the residual for the continuity equation. Similar results 
are obtained for the momentum and energy equations. The refore, the residuals on a fine grid can be transferred 
to a coarser grid so that the integral of the fluxes over the cell boundary is conserved simply by summing the 
fine-grid residuals that make up the coarse grid. 


(C12) 
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